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v^ ^ Abstract. We extend a template-based approach for synthesizing switch- 

ing controllers for semi-algebraic hybrid systems, in which all expres- 
k^j ' sions are polynomials. This is achieved by combining a QE (quantifier 

'^ elimination)-based method for generating continuous invariants with a 

qualitative approach for predefining templates. Our synthesis method 
CZ2 ' is relatively complete with regard to a given family of predefined tem- 

^ plates. Using qualitative analysis, we discuss heuristics to reduce the 

numbers of parameters appearing in the templates. To avoid too much 

human interaction in choosing templates as well as the high computa- 

^ ■ tional complexity caused by QE, we further investigate applications of 

lO ' the SOS (sum-of-squares) relaxation approach and the template poly- 

Cn , hedra approach in continuous invariant generation, which are both well 

y5^ ■ supported by efficient numerical solvers. 

^' 

O ■ 1 Introduction 

m 



Hybrid systems, in which computations proceed both by continuous evolutions 
as well as discrete jumps simulating transition from one mode to another mode, 
are often used to model devices controlled by computers in many application do- 
^ , mains [I] . Combining ideas from state machines in computer science and control 

Cy_' theory, formal analysis, verification and synthesis of hybrid systems have been 

an important area of active research. In verification problems, a given hybrid 
system is required to satisfy a desired safety property e.g. that the temperature 
of a nuclear reactor will never go beyond a maximum threshold, as it may cause 
serious economic, human and/or environmental damage, thus implying that the 
system will never enter any unsafe state. A synthesis problem is harder given 
that the focus is on designing a controller that ensures the given system will 
satisfy a safety requirement, reach a given set of states, or meet an optimality 
criterion, or a desired combination of these requirements. 

Automata-theoretic and logical approaches have been primarily used for ver- 
ification and synthesis of hybrid systems |2I4I42| . In |4I42| , a general framework 
for controller synthesis based on hybrid automata was proposed, which relies on 



backward reachable set computation and fixed point iteration. Two restrictions 
of this approach are (i) the computation of backward reachable set is hard for 
most continuous dynamics, and (ii) termination of the fixpoint iteration pro- 
cedure cannot be guaranteed, even for those hybrid systems whose backward 
reachable sets are easily computable. Thus most of the research, e.g. [14], fo- 
cuses on overcoming the above two restrictions. 

Recently, a deductive approach for verification and synthesis based on con- 
straint solving was proposed in |11I27I26I38I20I37] . The central idea is to reduce 
verification and synthesis problems of hybrid systems to invariant generation 
problems, much like verification of programs. As proposed in |16ll5l3f] . if in- 
variants are hypothesized to be of certain shapes, then corresponding templates 
with associated parameters can be used and the invariant generation problem 
can be reduced to constraint solving over parameters by quantifier elimination. 
This methodology is used in [35] for synthesizing switching controllers meeting 
safety requirements, while in [40j, the approach is extended for satisfying both 
safety and reachability requirements. A common problem with template-based 
method is that it heavily relies on a user specifying the shape of invariants that 
are of interest, thus making it interactive and user driven, raising doubts about 
its scalability and automation. Besides, the inference rules for inductive invari- 
ants in |39|38|40] are sound and complete for classes of invariants, e.g. smooth, 
quadratic and convex invariants, but are not complete for generic semi-algebraic 
sets. 

Inspired by |17|4|38j and [21], we extend in this paper the template-based 
invariant generation approach for synthesizing switching controllers of hybrid 
systems to meet given safety requirements. The paper makes the following con- 
tributions: 

— We formalize the solution to switching controller synthesis problem (Problem 
[T]) of hybrid systems in terms of continuous invariants (see Theorem [T]), 
and thus lay the foundation of the synthesis method based on continuous 
invariant generation using templates and constraint solving. 

— In the QE-based synthesis framework we use the method for continuous 
invariant generation proposed in j21| as an integral component. This method 
is proved in [21] to be sound and relatively complete with respect to a given 
shape of invariants (i.e. a given family of predefined templates). As a result, 
in contrast to the methods used in [39|38I40] , there is more flexibility in our 
approach because of the possibility of discovering all possible invariants of 
the given shape. 

— Using the qualitative approach proposed in [17] for analyzing continuous 
evolution in certain modes of a hybrid system, we can develop heuristics to 
determine a more precise shape of templates to be used as invariants, thus 
reducing the numbers of parameters appearing in templates. 

— We further improve the degree of automation and scalability of the template- 
based method in two ways: (i) for general polynomial templates, using sum- 
of-squares (SOS) relaxation, the constraint on parameters appearing in tem- 
plates is transformed into a semi-definite program{SDP), which is convex and 



thus can be solved efficiently; (ii) for linear systems and a special type of 
templates — template polyhedra, again by sacrificing relative completeness, 
the continuous invariant generation problem can be reduced to a BMI (bilin- 
ear matrix inequality) feasibility problem, which is also much easier to solve 
(numerically) than QE. 

Related Work 

Our work in this paper resembles [35] but differs in that: i) our method is cast in 
the setting of hybrid automata and searches for a family of continuous invariants 
that refine the original domains, rather than a single global controlled invariant; 
ii) a sound and complete criterion is used in continuous invariant generation; iii) 
various techniques are applied for scalability. 

The SOS relaxation approach has been successfully used in safety verifica- 
tion of hybrid systems. In |28l29j . the authors used the SOSTOOLS software 
package [30) to compute barrier certificates for polynomial hybrid systems. In 
|19I45| . the authors proposed a hybrid symbolic-numeric approach to compute 
exact inequality invariants of hybrid systems, by first solving (bilinear) SOS pro- 
gramming numerically and then applying rational vector recovery techniques. 

A necessary and sufficient condition for positive invariance of convex polyhe- 
dra for linear continuous systems was provided in [7] . This condition is extended 
to linear systems with open polyhedral domain for our need in the paper. Tem- 
plate polyhedra was used in |33I32) to compute positive invariants of hybrid 
systems by policy iteration, which differs from our treatment of the problem 
using BMI. Recently, a method for computing polytopic invariants for polyno- 
mial dynamical systems using template polyhedra and linear programming was 
proposed [55] . 

Mathematical programming techniques and relevant numerical solvers have 
also been widely applied to static program analysis. Actually, the template poly- 
hedra abstract domain was first proposed in [33] to generate linear program in- 
variants using linear programming. In [5], to verify invariance and termination 
of semi-algebraic programs, verification conditions are abstracted into numeri- 
cal constraints using Lagrangian relaxation or SOS relaxation, which are then 
resolved by efficient SDP solvers. 

In our recent work [IS] , we studied an optimal switching controller synthesis 
problem arising from an industrial oil pump system with piece-wise constant 
continuous dynamics. A hybrid approach combining symbolic computation with 
numerical computation was developed to synthesize safe controllers with better 
optimal values. 

Paper Structure. To be completed. 
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2 Problem Description 

Following |4I42| . we use hybrid automata to model hybrid systems. 

Definition 1 (Hybrid Automaton). A hybrid automaton (HA) is a system 
n = {Q,X,f,D,E,G), where 

• Q = [qi, . . . , g„j} is a set of discrete states; 

• X = {xi, . . . , Xn} is a set of continuous state variables, with x ~ {xi, . . . ,Xn) 
ranging over M"; 

• f '■ Q ^ (M" — ?> M") assigns to to each discrete state q Cz Q a vector field f^; 

• D : Q -^ 2^'' assigns to each discrete state q ^ Q a domain Dq C R"; 

• E <Z Q X Q is a set of discrete transitions; 

• G : E —>■ 2^ assigns to each transition e & E a switching guard Ge ^ K" • 

Remark 1. For ease of presentation, we make the following assumptions: 

— for all q ^ Q, iq is polynomial vector function; besides, fg defines a complete 
vector field, that is, for any xq G M", the solution to the differential x = fg 
uniquely exists on [0, oo); 

— for all g G Q and all e ^ E, Dq and Ge are closed semi- algebraic setqj; 

— the initial condition in each discrete mode is assumed to be identical with 
the domain, and all reset functions are assumed to be identity mappings. 

We use a nuclear reactor system discussed in [3112117] as a running example 
through this paper. 

Example 1. The nuclear reactor system consists of a reactor core and a cooling 
rod which is immersed into and removed out of the core periodically to keep the 
temperature of the core, denoted by x, in a certain range. Denote the fraction 
of the rod immersed into the reactor by p. Then the initial specification of this 
system can be represented using the hybrid automaton in Fig. [TJ 

The semantics of a hybrid automaton "H can be defined by the set of trajecto- 
ries it accepts. For the formal definitions of hybrid time set and hybrid trajectory 
the readers are referred to [32]- We denote the set of trajectories of % by Tr{'H), 
ranged over w, wi, . . .. All trajectories of T-l starting from the initial state (go, xq) 
is denoted by Tr{'H){qo,xo). 

The domain of a hybrid automaton H is defined as D-h = UgeQd'?} ^ Dq). 
We call 'H non-blocking if for any [q, x) g Du, there is a hybrid trajectory from 



^ A set A C R" is called semi-algebraic if there is a quantifier-free polynomial formula 
(p s.t. A = {x e R" 1 (^(x) is true} . 



qi: no rod 




G41 



p = 



x^x/10-6p-50 
04 = 0<p<l 



G34. 



q2: being immersed 

X = a;/10-6p-50 
D2 = <P <1 




g4 : being removed 53 : immersed 

Fig. 1. Nuclear reactor temperature control. 



(g, x) which can either be extended to infinite time t = cxo or execute infinitely 
many discrete transitions; otherwise % is called blocking. 

A safety requirement S assigns to each mode q ^ Q a safe region Sq C M", i.e. 
S = Ugegd?} ^ ^q)- Alternatively, there could be a global safety requirement 
S which all modes are required to satisfy. 

According to [1] , the switching controller synthesis problem with regard to a 
given safety requirement can be formally defined as follows: 

Problem 1 (Controller Synthesis for Safety). Given a hybrid automaton % and 
a safety property 5, find a hybrid automaton Ti.' = (Q, X, f, D', E, G") such that 

(rl) Refinement: for any q (z Q, D' C Dq, and for any e ^ E, G'^ C Ge', 

(r2) Safety: for any trajectory u that T-L' accepts, if ((7,x) is on uj, then x G Sq] 

(r3) Non-blocking: H' is non-blocking. 



If such n' exists, then SC = {G^. C I 
satisfying safety requirement S, and D-w 
invariant set rendered by SC. 



I e G E} is a switching controller 
[JqGQiil} ^ ^q) is th^ controlled 



3 A QE-Based Approach 



3.1 Continuous Invariant 

Along the line of f38^ , we consider the switching controller synthesis problem by 
combining a relative complete method for generating continuous invariants in 
pi] , and heuristics for predefining templates for these invariants using qualitative 
analysis in |17) . Below, we review the concept of continuous invariant used in 
[2T] based on a related concept in [26] . 

Definition 2 (Continuous Invariant (CI)). Given a mode q ^ Q in a hybrid 
automaton %, a set P C M" is called a continuous invariant of (Dq,fq) if for 



all xq G P n Dq and all T > 0, the solution x(ijj of x. — fg(x) over [0,T] with 
x(0) = xq satisfies 



(Vte[o,r].x(i)ez?, 



(Vie [0,T].x(t) eP) 



Intuitively, P is a CI of {Dq, fg) if any continuous evolution starting from the 
intersection of P and Dq stays in P as long as it is still in Dq. If Dq ~ K", then 
a CI of {Dq,{q) coincides with the standard (positive) invariant set (see [S]) of 
the dynamical system defined by f^; otherwise if Dq is a proper subset of R", 
then generally the notion of CI is weaker. 

Example 2. Suppose Dq = a: > and iq = {—y, x). Obviously P = y > is not 
a positive invariant set of fg, whereas P is a CI of {Dq, iq) according to Definition 
m See Fig. [5] for an illustration. 
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Fig. 2. Illustration of continuous invariant. 



3.2 The Abstract Synthesis Procedure 

To solve Problem [T] amounts to refining domains and guards of H. by removing 
had states from domain D-^. A state ((7,x) S D-^ is 6ad if the hybrid trajectory 
starting from {q, x) either blocks 'K or violates 5; otherwise, it is called a good 
state. From Definition [2l we observe that the set of good states of V, can be 
approximated using continuous invariants, which results in the following solution 
to Problem [1] 

Theorem 1. Let T-i and S be the same as in Probleml^ Suppose D' is a closed 
subset o/K" for all q £ Q and D' ^ for at least one q. If we have 

(el) for all qeQ, D' CDqH Sq; 



We assume the existence and uniqueness of solutions is guaranteed by f. 



(c2) for all q E Q, D' is a continuous invariant of {Hq,iq) with 

H, = { u G0^ 

e={q,q')eE 

where G'^ =^ GeC] D' , and A'^ denotes the complement of A in R", then the HA 
n' = {Q,X,f,D',E,G') is a solution to ProblemUi 

Proof. Please refer to the Appendix [Xj D 

Remark 2. Intuitively, by (cl), D' is a refinement of Dq and is also contained in 
the safe region; by (c2), any trajectory starting from D' will either stay in D' 
forever, or finally intersect one of the transition guards enabling jumps from q 
to a certain g', thus guaranteeing satisfaction of the non-blocking requirement. 

Based on Theorem [U we give below the steps of a template-based method 
for synthesizing a switching controller. 

(si) Template assignment: assign to each q G Q a template parametrically 
specifying D' which can be seen as a refinement of Dq and will be instan- 
tiated to be the continuous invariant at q: 

(s2) Guard refinement: refine guard Gg for each e = {q,q') E E hy setting 
Gg = Ge n Uq, ; 

(s3) Deriving syntiiesis conditions: encode (cl) and (c2) in Theorem [T] into 
constraints on parameters appearing in templates; 

(s4) Constraint solving: solve the constraints derived from (s3) in terms of the 
parameters; 

(s5) Parameters instantiation: find an appropriate instantiation of D' and 
Gg such that D' are closed sets for all q G Q, and D' is nonempty for at 
least one g e Q; if such an instantiation is not found, we choose a new set 
of templates and go back to (si). 

Remarks: 

1. The implementability of the above method depends on the language used 
to specify the hybrid system, the safety property, as well as the templates 
chosen for their refinements. If all appearing expressions are specified using 
polynomials, the computability of the abstract procedure is guaranteed by 
Tarski's result [41 . This will be assumed in the rest of this paper. 

2. In (s3), condition (cl) can be encoded into a first-order polynomial formula 
straightforwardly; encoding of (c2) into first-order polynomial constraints is 
based on our previous work in |21) about a relatively complete method for 
generating CIs (see Section 15^ . 

3. We use quantifier elimination (QE) to solve the first-order polynomial con- 
straints obtained in (s4). 

4. The shape of chosen templates in (si) determines the likelihood of success of 
the above procedure, as well as the complexity of QE in (s4). In Section [BTH 
we discuss heuristics for choosing appropriate templates using the qualitative 
analysis discussed in |17] . 



3.3 A Relatively Complete Method for Generating CIs 



In |21) we presented a sound and relatively complete approach for generating 
semi-algebraic CIs for (Dq^tq) with semi-algebraic Dq and polynomial vector 
function f^. We review here the key ideas; for details the reader can consult [3T]. 
Below, we drop the subscript corresponding to the mode q. 

The basic idea can be explained as follows for the simplest case, namely 
D = /i(x) > and P = p(x) > 0. Let dP = p(x) = be the boundary of P and 
x(i) be the continuous evolution of f starting from xq. It can be shown that P 
is a CI of (D, f ) if and only if for all xq £ SF n Z? : 



3e>0VtG [0,e].p(x(t)) eP 



(1) 



Intuitively, ([T]) means that trajectories starting at the the boundary of P will 
stay in P for a small amount of time. 

Given that p and f are polynomials and thus analytic, the Taylor expansion 
ofp(x(t)) at t = 



dp dp t 

^W^))-^(^")+dI-^+dZ^-2!^ 

converges in a neighborhood of 0. 

Define the Lie derivatives of p along f , L^p : R" — > 

• ifP(x) = p(x) ; 

• ^pix) = (|jLrV(x),f(x)), for n > 0, 

where (■, •) is the inner product of two vectors, i.e. ( (ai, 
X]r=i o-i^i- Using Lie derivatives, ([2|) rewritten as 



(2) 



for n G N, as follows: 



p(x(i)) = L?p(xo) + L}p{xo) ■ t + LIp{ko) 



2! 



, On), {h 






,bn)] 



(3) 



Combining ^ and ([3]), our main result of continuous invariant generation 
(in the simplest case) can be stated as follows. 

Theorem 2 (Necessary and Sufficient Criterion for CIs |21| ). Given a 
system {D,{) with D = /i(x) > 0, it has a continuous invariant of the form 
P = p(x) > if and only if Vx.(p(x) = A h{x.) > — > ip{p, f)), where 



^(p,f) 



L}p{x) > 

V ifp(x) = A ifp(x) > 

V ■•• 

V L^p(x) = A • • ■ A Lf"'" V(x) = A if^'Vlx) > 

V Ljp(x) = A ■ • • A if"'" V(x) = A Lf^'pix.) = 



with Np_f G N computed from p and f . 
Proof. Please refer to |21) . 



Intuitively, Theorem [5] means that on the boundary of P, up to the A'p^f-th 
order, the first non-zero higher order Lie derivative of p w.r.t f is non-negative. 

The above theorem can be generahzed for parametric polynomials p(u,x), 
thus enabling us to use polynomial templates and QE to automatically discover 
CIs. Such a method for CI generation is relatively complete, that is, if there 
exists a CI in the form of the predefined template, then we are able to find one. 

Example 3. Consider the system (]R^, f) from [321 with f = (i = 1 — y, y = x), 
which has a continuous invariant p > with p=— (—x^ ^ 2/^ + 2y)^, defining 
the circumference of a circle. 

In [39 , sound and complete inference rules are given for invariants that are 
linear, quadratic, smooth or convex. However, it was pointed out in |39j that all 
these rules failed to prove the invariance property of p > 0, as p is not linear 
or quadratic, nor is it smooth or convex. Furthermore, by a simple computation 
we get L^p = for all fc > 1, so the sound but incomplete rule in 



which involves only strict inequalities over finite-order Lie derivatives is also 
inapplicable. However, from L\p = we get Npf = 0, and then according to 
Theorem[21 p > can be verified since VxVy. (— (— x^ — 2/^ + 2y)^ = — > true) 
holds trivially. 

Although the rule in [26] can also be used to check the invariant p > 0, 
generally it only works on very restricted invariants. Even for linear systems like 
(R, a; = x), it cannot prove the invariant x > because Vx.x > is obviously 
false, while our approach requires \/x.{x = — >■ true) which is trivially true. 

The above examples show the generality and flexibility of our approach, using 
which it is possible to generate CIs in many general cases, and hence gives more 
possibility to synthesize a controller based on our understandings of the kind of 
controllers that can be synthesized using methods in [38;40p37j . 

3.4 Heuristics for Predefining Templates 

The key steps of the qualitative analysis used in IXTl are as follows. 

1. The evolution behavior (increasing or decreasing) of continuous variables in 
each mode is inferred from the differential equations (using flrst or second 
order derivatives); 

2. control critical modes, at which the maximal (or minimal) value of a contin- 
uous variable is achieved, can be identified; 

3. the safety requirement is imposed to obtain constraints on guards of transi- 
tions leading to control critical modes, and 

4. then this information is propagated to other modes. 

Next, we illustrate how such an analysis helps in predefining templates for 
the running example. 

Example 4 (Nuclear Reactor Temperature Control). Our goal is to synthesize a 
switching controller for the system in Example [T] with the global safety require- 
ment that the temperature of the core lies between 510 and 550, i.e. Si = 510 < 
a: < 550 fori = 1,2, 3, 4. 



1) Refine domains. Using the safety requirement, domains Di for z — 1, 2, 3, 4 
are refined by Df ^ A n S^, e.g. Z?^ = p = A 510 < a; < 550 . 

2) Infer continuous evolutions. Let li = x/10 — 6p — 50 = be the zero-level 
set of X and check how x and p evolve in each mode. For example, in Z?2j 
i > on the left of li and i < on the right; since p increases from to 1, a; 
first increases then decreases and achieves maximal value when crossing li. 

3) Identify critical control modes. By 2), q2 and g4 are critical control 
modes at which the evolution direction of x changes. 

4) Generate control points. By 3), we can get a control point i?(5/6, 550) 
at q2 by taking the intersection of h and the safety upper bound x — 550; 
and F(l/6, 510) can be obtained similarly at (74. 

5) Propagate control points. E is backward propagated to A{0, a) using 
trajectory ae at ^2: and then to C(l, c) using trajectory CA at 94; similarly, 
by propagating F we get D and B. (See Fig. |31) 

6) Construct templates. For brevity, we only show how to construct Dj- 
Intuitively, p — 0^ p — 1, ae and bd forms the boundaries of D'2- In order to 
get a semi-algebraic template, we need to fit ae and bd by polynomials using 
points A, E and B, D respectively. By 2), ae has only one extreme point E 
in Z?| and is tangential to a; = 550 at E. The simplest algebraic curve that 
can exhibit a shape similar to ae is the parabola through A, E opening 
downward with h = P = ^ the axis of symmetry. Therefore to minimize the 
degree of terms appearing in templates, we do not resort to polynomials with 
degree greater than 3. This parabola can be computed using the coordinates 
of ^,^ as: a;-550-f|(a-550)(p-|)2 = 0. 





0^ 0.4 0.6 



Fig. 3. Control points propagation. 



Through the above analysis, we generate the following templates: 

• D[^ p = A r3lO<x<a; 

• D'2^ 0<p<l A x-b>p{d-b) A x~550-§ia- 550)(p - 1)^ < ; 

• D^ = p = 1 A d < x < 550 ; 

• D'^= 0<p<lAx-a<p{c-a) A a; - 510 - ||(«i - 510)(p - i)^ > , 
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in which a, 6, c, d are parameters. Note that without qualitative analysis, a single 
generic quadratic polynomial over p and x would require { t ) = ^ parameters. 

The above heuristics works well on planar systems and can also be applied 
to three-dimensional systems. We are further generalizing the heuristics to cover 
a wider class of hybrid automata. 

Based on the framework presented in Section 13.21 we show below how to 
synthesize a switching controller for the system in Example 2] step by step. 

Example 5 (Nuclear Reactor Temperature Control Contd.). 



(si) The four invariant templates are defined in Section [3.41 
(s2) The four guards are refined by setting G[j = Gij Ci D'j : 

• G'i2 = p = 0A6<a;<a; 

• G"23 = p = 1 A d < X < 550 ; 

• ^34 = p=lAd<x<c; 

• G^i ^ p = A 510<x< a. 

(s3) Using D'^ and G^ ■ we can derive the synthesis condition, which is a first- 
order polynomial formula in the form of = VxVp.(^(a, 6, c, d, a;,p). We do 
not include (j) here due to its big size. 

(s4) By applying QE to <j) we get the following solution to the parameters: 



6575 
12 



A b 



4135 



A c 



4345 



A d 



6145 
12 



(4) 



(s5) Instantiate D^ and G'^j by (g]). It is obvious that all D- are nonempty closed 
sets. According to Theorem [TJ we get a safe switching controller for the 
nuclear reactor system. The left picture in Fig. |4]is an illustration of Di^. 



Fig. 4. Shape of synthesized continuous invariants. 



In [T7], an upper bound x = 547.97 for G12 and a lower bound x — 512.03 
for G34 are obtained by solving the differential equations at mode (72 and 54 
respectively. By Q, the corresponding bounds generated here are x < 
547.92 and X > ^ = 512.08. 



6575 
12 
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As should be evident from the above resuhs, in contrast to 17 , where differ- 
ential equations are solved to get closed-form solutions, we are able to get good 
approximate results without requiring closed-form solutions. This indicates that 
our approach should work well for hybrid automata where differential equations 
for modes need not have closed form solutions. 



4 Synthesis by Generating CIs Numerically 

The QE-based approach crucially depends upon quantifier elimination tech- 
niques. It is well known that the complexity of a general purpose QE method 
over the full theory of real-closed fields is doubly exponential in the number 
of variables [5]. Therefore it is desirable to develop heuristics to do QE more 
efficiently. As shown in Section 13.41 qualitative analysis helps in reducing the 
number of parameters in templates. Another possible way to address the issue 
of high computational cost is resorting to numerical methods. In this section, 
we will discuss the application of two such approaches to the nuclear reactor 
example. 

4.1 The SOS Relaxation Approach 

Let IR[a;i, X2, ■ ■ ■ , a;„], or M[x] for short, denote the polynomial ring over variables 
Xi,X2t ■ ■ ,Xn with real coefficients. A monomial is a special polynomial in the 
form of Xi^X2^ ■ ■ ■ x^" with {ai,a2, ■ . ■ , a„) G N". Any polynomial p(x) e R[x] 
of degree d can be written as a linear combination of ("^ ) monomials, i.e. 



PV^J ~ / , C(q,j Q,2,...,a„) ' • 



ai+a2H l-an<d 



A polynomial p is called an SOS (sum-of-squares) if there exist s polynomials 
qi,q2,...,qs s.t. 

p- E'^^ 

l<i<s 

It is obvious that any SOS p is non- negative, i.e. Vx e M".p(x) > . 

The basic idea of SOS relaxation is as follows: to prove that a polynomial p is 
nonnegative, we can try to show that p can be decomposed into a sum of squares, 
a trivially sufficient condition for non- negativity (but generally not necessary); 
similarly, to prove p > on the semi-algebraic set g > 0, it is sufficient to find 
two SOS ri, r2 such that p = ri + r2 • q. 

SOS relaxation is attractive because the searching for SOS decomposition 
can be reduced to a semi-definite programming (SDP) problem according to the 
following equivalence |25) : 

A polynomial p of degree 2d is an SOS if and only if there exists a semi- 
definite matrix Q such that p = q-Q-q^, where q is a ("^ j-dimensional 
row vector of monomials with degree < d. 
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SDP is a convex programming that is solvable in polynomial time using numerical 
methods such as the interior point method [H]. Therefore the SOS computation 
is a tractable problem. 

We now show how SOS can be related to CI generation. Let p > be a 
parametric template defined for the system {h > 0, f ). By Theorem[2l a sufficient 
condition for p > to be a CI of {h > 0, f ) is 

Vx.(p(x) = A h{x.) > — > LIp{x) > 0) , 

which can be further strengthened to 

Vx. (/i(x) > ^ LIp{x) > 0) . (5) 

Using SOS relaxation, a sufficient condition for ([5]) is 

LIp = Si + 32 ■ h + e , (6) 

where si,S2 are SOS and e is a positive constant. Solve (jG)) for the parameters 
appearing in p and then we can get a CI p > of (ft, > 0, f). 

For the nuclear reactor example, we define two general quartic templates 

p>OAp<lA J2 C(„i,a.)-P"'^"^<0 

and 

p>OAp<lA J2 ««(ai,a2)-p"'a^"' <0 

for mode 92 and 94 respectively. Using the SOS relaxation techniques discussed 
above, the following refined domains are obtained: 



D'2= 0<p<l A 34468.9- 9941. 89p+1114.38p2 + 67.3261p3 + 0.925p4 



• £>; = p = A 510 < a; < 547.85 ; 

129.316a; + 37.7294pa; - 4.9669^^2; - 0.1303^^2; + 0.1212a;2 - 0.0358px2 + 
0.0054p2j;2 < ; 
, D'^^ p^l A 512.09 < a; < 550 ; 

• D'4=0<p<lA 46082.8 + 8787.98p + 1473.0p2 + 93.8933^^ + 1.635^" - 
174.456a; - 34.1747pa; - 4.7397^^2; - 0.1829p3a; + 0.1649a;2 + 0.0332pa;2 + 
0.0037p2a;2 < . 

The picture in the middle of Fig. |4] illustrates the synthesized D'2. 

4.2 The Template Polyhedra Approach 

Polyhedral sets are a popular family of (positive) invariants of linear (continuous 
or discrete) systems 0. A convex polyhedron in R" can be represented using 
linear inequality constraints as Qx < p, where Q € M''^" is an r x n matrix, and 
jj.g]^nxi^^g]girxi g^j.g column vectors. 

Given a linear continuous dynamical system x — Ax with A e M"^", the 
following result about (positive) polyhedral invariant set is established in [7]. 
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Proposition 1. The polyhedron Qx < p is a positive invariant set ofx. = Ax if 
and only if there exists an essentially non-negativ^ matrix H G W^"^ satisfying 
HQ = QA and Hp < 0. 

By simply applying the famous Farkas ' lemma jll) , we can generalize Propo- 
sition [T] and give a sufficient condition for polyhedral CIs of linear dynamics with 
open polyhedral domain. 

Proposition 2. Let f S Ax + b and D = ex < a, where a G M, b e R"^^ is a 
column vector, and c G K^^" is a row vector. Then the polyhedron Qx < p is a 
CI of the system {D, f ) if there exist essentially non-negative matrix H G R^^*^, 
and non-negative column vectors ?? > 0, ^ > 0, A > inW^^ such that 

(1) i7Q + ec-diag(A)g.4 = ; 

(2) Hp + r] + £_a + diag(A)Qb = ; 

(3) ( + V>0 , 

where diag(A) denotes the r x r diagonal matrix with the main diagonal A. 

Proof. Please refer to Appendix |B] D 

Proposition[2]serves as the basis of automatic generation of polyhedral CIs for 
linear systems. To reduce the number of parameters in a polyhedral template, 
we propose the use of template polyhedra. The idea is to partly fix the shape 
of the invariant polyhedra by fixing the orientation of their facets. Formally, a 
template polyhedron is of the form Qx < p where Q is fixed a priori and p is to be 
determined. Any instantiation of p from W^^ produces a concrete polyhedron. 
In this paper, since the system is planar, we choose Q in such a way that its row 
vectors form a set of uniformly distributed directions on a unit circle, i.e. 

i ~ \ i — 1 
Qi — (cos( 27r),sin( 27r)) 

for 1 < i < r, where q^ denotes the i-th row of Q. It is easy to see that Qx < p 
is just a rectangle when r = 4, and an octagon when r = 8. 

In order to determine p, we have to solve the constraints derived from Propo- 
sition [2j Note that since both H and p are indeterminate, the constraint (2) 
becomes bilinear, making the problem NP-hard ^43j to solve. It is however still 
more tractable using modern BMI (bilinear matrix inequality) solvers compared 
to QE. The details of applying numerical solvers will be discussed in the Con- 
clusion part. 

Using octagonal templateijj for mode (72 and (74 in the nuclear reactor exam- 
ple, we obtain the following refined domains. 



^ A square matrix is essentially non-negative if all its entries are non-negative except 
for those on the diagonal. Besides, given a matrix M, in this paper the notation 
M > 0, M > and M = should be interpreted entry-wisely. 

** To reduce the number of facets needed in the template, we scaled the variable x by 
a factor of 0.2, i.e. let x = 5x' , and rescaled the generated invariants by 5. 
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• D[ = p = A 510<x < 545.50 ; 

• D'2= Q{p, xY < Pi with 



. 0.0000 0.7071 1.0000 0.7071 0.0000 -0.7071 -1.0000 -0.7071 

and 

pi = (5.0000 392.4429 549.9276 385.8688 0.0000 -367.6169 -514.5244 -360.2745)'^; 

• D^ = p = 1 A 514.50 < a; < 550 ; 

• D'^ = Q{p, x)'^ < p2 with the Q in D2 and 

P2 — (5.0000 384.1267 545.5548 384.7484 0.0000 -360.6299 -510.1431 -360.1948)"^. 

In Fig. m the picture on the right illustrates the synthesized D2. 



5 Conclusion and Discussion 

We have extended a template-based approach for synthesizing switching con- 
trollers for semi-algebraic hybrid systems by combining symbolic invariant gen- 
eration methods using quantifier elimination with qualitative methods to deter- 
mine the likely shape of invariants. We have also investigated the application of 
numerical methods to gain high level of scalability and automation. A summary 
comparison of the three proposed approaches, i.e. the QE-based, SOS-relaxation 
and template-polyhedra approaches, can be given in the following aspects. 

— Applicability: the QE-based approach can be applied to any semi-algebraic 
system; SOS relaxation techniques can be applied to semi-algebraic systems 
for which SOS encoding is possible; the template-polyhedra approach is only 
applicable to linear systems. 

— Design of Templates: the QE-based approach demands much heuristics 
in determining templates, while the other two need little human effort. 

— Relative Completeness: only the QE-based approach is relative complete 
w.r.t. the predefined family of templates, but we believe that the template- 
polyhedra approach can be made relatively complete by improving Prop. [H 

— Quality of Controllers: the QE-based and SOS-relaxation approaches can 
generate arbitrary (non-convex) semi-algebraic invariants, while the template- 
polyhedra approach can only generate convex polyhedral invariants; for the 
nuclear reactor example, we can see from Fig.|3]that the QE-based approach 
produced larger refined domains and transition guards, but such superiority 
is not a necessity and relies greatly on the quality of heuristics. 

— Computational Cost: for the QE-based approach, we have used the alge- 
braic tools Redlog Hm and QEPCAD (the sifq function) [5] to perform QE; 
for the numerical approaches, we use the MATLAB optimization toolbox 
YALMIP |23l24j as a high-level modeling environment and the interfaced 
external solvers SeDuMi [36] and PENBMI [lS\ (the TOMLAB [13] version) 
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Table 1. Templates and time cost of three controller synthesis approaches. 



Approach 


QE-based 


SOS-relaxation 


template-polyhedra 


Tool 


Redlog + sifq 


YALMIP + SeDuMi 


YALMIP + PENBMI 


Template 


NR 


quadratic, #PARMS — 4 


generic quartic 


8 facets 


TS 


quadratic, #PARMS — 2 


generic quartic 


10/12 facets 


Time 

(sec) 


NR 


12.663 


1.969 


0.578 


TS 


7.092 


1.609 


1.437 



to solve the underlying SDP and BMI problems respectively. Table [U shows 
the time cost of three approaches applied to the nuclear reactor (NR) ex- 
ample as well as a thermostat (TS) example from [14]. All computations 
are done on a desktop with a 2.66 GHz CPU and 4 GB memory. We can 
see that for these two examples the QE-based approach is consistently more 
expensive in time compared to numerical approaches. 
— Soundness: the QE-based approach is exact while the other two approaches 
suffer from numerical errors which would cause the synthesis of unsafe con- 
trollers. The justification for use of numerical methods is that verification 
is much easier than synthesis. For example, we have verified posteriorly and 
symbolically the controllers synthesized by both numerical approaches in 
this paper. We could also directly encode some tolerance of numerical errors 
into the synthesis constraints to increase robustness and reduce the risk of 
synthesizing bad controllers. 

Our analysis of a nuclear reactor example suggests the effectiveness of all 
three proposed approaches. We are currently experimenting with these (and 
more other) methods on more complex examples. We believe that there exists 
no single method that can solve all the problems. A practical way is to select 
the most suitable one(s) for any specific problem. 

Although the focus of this paper is on the switching controller synthesis 
problem subject to safety requirements, we plan to extend the proposed approach 
for reachability and/or optimality requirements as well, by incorporating our 
previous results on asymptotic stability analysis |22j and a case study in optimal 
control l46l. 
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A Proof of Theorem [T] 

We need the following definitions [55] to prove Theorem [T] 

Definition 3 (Hybrid Time Set). A hybrid time set is a sequence of intervals 
T = {Ii}iL() (^ can be cx)j such that: 

• li ~ [Ti, Tj'] with Ti < t[ ~ T,;+i for all i < N ; 

• if N < oo, then I^ — [tn,t'j^) is a right-closed or right-open nonempty 
interval (V^ may be cxdJ; 

. To = . 

Given a hybrid time set, let (t) = N and ||t|| = ^,-o('''i' ^ '''i) ■ 

Definition 4 (Hybrid Trajectory). A hybrid trajectory of TL starting from 
an initial point (qo,Xo) G Df{ is a triple uj = (t, a,/3), where r = {/i}^o *■' '^ 
hybrid time set, and a — {ai : li -^ Q}iLo and /3 = {Pi : li — ?> M"}flg are two 
sequences of functions satisfying: 

1. Initial condition: ao[0] = go o,nd /3o[0] — xq; 

2. Discrete transition: for alii < (r), e = (^^(t,'), ai+i(Ti+i)) G E, (ii{T[) G Ge 
and /3i+i(Ti+i) = /3j(t/); 

5. Continuous evolution: for all i < (r) with Ti <t[, if q — ai{Ti), then 

(1) for all t e h, ai{t) = q, 

(2) Pi (t) is the solution to the differential equation x = f ^ (x) over li starting 
from Pi{Ti), and 

(3) for all t€ [T,;,T/),/3,(i) & Dq . 

Proof of Theorem [l] 

Proof. We prove that the three requirements in Problem [T] are satisfied by %' . 

(rl) By (cl), we get D'^ C Dq for aU g G Q; by the definition of G'^, G^ C Ge 
for all e e E. 

(r2) Suppose lu — (T,a,P) is a hybrid trajectory starting from (gojXo) G -D-h'- 
We prove 

Vz<(r}V^G/,.A(i)G5„,(t) (7) 

by induction on (r). 

If (t) = 0, then Iq = [0,Tq) is a right-open or right-closed interval for some 
To > 0. If To = then /q = {0}. By (cl) and condition 1 of Definition H we have 
/3o(0) = xo G I?^p C Sqg — Sao^Q). If To > 0, by condition 1 and 3 of Definition 
Has weh as (cl), we have for aU t G [0,To), /3o(i) G D'^^^ C Sq,, == S'„o(t)- If 
Jo ^ [0,To), i.e. Iq = [0,To], by noticing that Dq^ is a closed set and /3o(i) is 
continuous over Iq, we get Po{To) G D'^^ C Sq^ = 'S'qo(To)- Thus we have proved 
in all cases, Vi G Io-Po{t) G ^'^^(f). 

Assume ([7]) holds for (r) = fc > 0. When (r) = fc + 1, by assumption, for 
all i < A; and all t e li we have Pi{t) G Sa-(t)- By condition 2 in Definition 
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El there exists e = {qk,qk+i) e E such that afc(T(.) = g^, afc+i(rfc+i) = qk+i 
and /3fc+i(Tfe+i) = /3fc(T^) £ Gg C D'^^^_^. Consider the hybrid trajectory starting 
from ((7fe+i,/3fe+i(T/c+i)) G i'w- By applying the same analysis we do for case 
(r) = 0, we can get /3fc+i(i) G S^^^^i^t) for all i G Ik+i- Thus we have proved 
that O holds for (r) = fc + 1. Then by induction, © holds for all (r). 

(r3) Given (go,xo) G -D-h' (hence D' ^ 0), we will construct a non-blocking 
hybrid trajectory starting from (qo,xo). 

Suppose x(t) is the continuous evolution defined by fg„ starting from Xq. Let 
Tmax be the maximal positive T satisfying 

x(i) G D'^^ n Hq, for all t G [0, T) , (8) 

if such T exists, and Tmax = otherwise. 

If Tlnax = oo, then we already get an infinite hybrid trajectory starting from 
(90, xo). 

If Tlnax < OO, then by the completeness of f^^, x(i) must exist on [0, T^aax + s] 
for some e > 0. We assert that 

x(r^ax)G (il,„r- U G', . (9) 

e={qo,q')<^E 

If not, i.e. x(rmax) G Hq^, then there exists < e' < £ s.t. x(i) G Hq^ on 
[0,Tniax + s'Ij because by assumption Hq^ is an open set. Then by (c2) and the 
definition of continuous invariant, we get x(t) G D' n Hq^ on [0,Tmax + £']: so 
T'max could not bc maximal. Therefore (jH]) holds. Then there exists e — {qq, q') G 
E such that x(Tiiiax) G Gg C D', , so we can make a discrete jump from qo to q' 
and extend the hybrid trajectory by continuous evolution at q' . 

Such extension either ends with a trajectory satisfying |[t|| = cxd or goes on 
forever resulting (t) — oo. D 



B Proof of Proposition [2] 

To prove Proposition [51 we first introduce the famous Farkas' lemma [11] in the 
theory of linear programming. 

Lemma 1 (Farkas' Lemma). For linear formulas Pj,rk G K[x], the formula 
/\iejPj > ^ AfceK^fe > is unsatisfiahle (over the reals) if and only if there 
exist non-negative constants fi, fij {j G J) and Vk {k G K) such that 

jeJ keK 

and at least one of fij,fj, is strictly positive. 
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Proof of Proposition [2] 

Proof. Let q^ denote the i-th row vector of the matrix Q and pi denote the 
i-th entry of column vector p. Then q^x = pi stands for the i-th facet of the 
polyhedron Qx < p. By a generalization of Theorem [5] (see [H] for the detail), 
Qx < p is a CI of (ex < a, Ax + b) if for all 1 < i < r, the following implication 
holds: 

Qx < p A q^x = Pi A ex < a =^ <li{Ax + b) < , 
which is equivalent to 

- Qx + p > A q^x - Pj > A qj(Ax + b) > A -ex + a > (10) 

is unsatisfiable. 

By Lemma [Tl the unsatisfiability of p^ is equivalent to the existence of 
constant 7*, and non-negative constants 7j\7f , . . . ,7|~\7]+\ . . . 7[, 77^, Aj, S,i 
such that 

Y^ j,{-qiX + Pi) + Xi(li{Ax + h)+i,{-cx + a) + T]^ = (11) 

l<i<r 

and £,i + rji > 0. By equating the coefficients of the left side of ([TT|) to 0, we get 

(1) ELi(-7.qO + A,q,^-e,c = 0; 

(2) Xir=i 7*Pi + Ajq^b + da + 77^ = ; and 

(3) d+m>0, 

the matrix form of which corresponds to the three conditions in Proposition [51 

D 
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